Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 3 November 2010 (MN M£X style file v2.2) 



The formation of disc galaxies in a ACDM universe 
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ABSTRACT 

We study the formation of disc galaxies in a fully cosmological framework using adap- 
tive mesh refinement simulations. We perform an extensive parameter study of the 
main subgrid processes that control how gas is converted into stars and the coupled 
effect of supernovae feedback. We argue that previous attempts to form disc galaxies 
have been unsuccessful because of the universal adoption of strong feedback combined 
with high star formation efficiencies. Unless extreme amounts of energy are injected 
into the interstellar medium during supernovae events, these star formation parame- 
ters result in bulge-dominated SO/Sa galaxies as star formation is too efficient at z ~ 3. 
We show that a low efficiency of star formation more closely models the sub-parsec 
physical processes, especially at high redshift. We highlight the successful formation 
of extended disc galaxies with scale lengths r& — 4 — 5 kpc, flat rotation curves and 
bulge-to-disc ratios of B/D~ 1/4. Not only do we resolve the formation of a Milky 
Way-like spiral galaxy, we also observe the secular evolution of the disc as it forms 
a pseudo-bulge. The disc properties agree well with observations and are compatible 
with the photometric and baryonic Tully-Fisher relations, the Esfr - S gas (Kennicutt- 
Schmidt) relation and the observed angular momentum content of spiral galaxies. We 
conclude that the underlying small-scale star formation physics plays a greater role 
than previously considered in simulations of galaxy formation. 

Key words: galaxies:evolution - galaxies:formation - galaxies: haloes - galaxies:spirals 



1 INTRODUCTION 

The prevailing pic ture of galaxy forma ti on emerged more 
than 30 yr ago (|White fc Reesl 1 19781 ; iFall fc Efstathioul 
1980i). Within the framework of the broadly accepted A Cold 
Dark Matter (ACDM) scenario (|Komatsu et alj|2009h . grav- 
ity assembles structures in a bottom-up fashion. Haloes of 
dark ma t ter a c quire angular momentum via tidal torques 
(|Peeblesl 1 19691 ; I Fall fc Efstathioul Il980h from interacting 
structures, and as gas cools and condenses into their cen- 
tral parts, star-forming galaxies form. A realistic angular 
momentum content can be accounted for if most of the an- 
gular momentum is retained in the assembly process. In 
this picture, the hos t halo is respon sible for the final galaxy 
characteristics (e.g. iMo et aljfl998| y While several aspects 
of the theory of galaxy formation are still being developed, 
e.g. the underlying physics of the miss ing satellite problem 
|Klvpin et al .111 9991 ; fMoore et alj|l999h and the role of cold 



Given the complexity and non-linearity of the in- 
volved processes, computer simulations have become the 
ideal tool for studying the formation of structure. The 
formation of a late-type spiral galaxy, such as our own 
Milky Way, has been studied numerically in fully ACDM 



cosmo log ical context by many a uthor s ( e.g. lAbadi et al 
| 2003bl ; ISommer-Larsen et al.l 120031; IGovernato et al 
120041; iRobertson et all 12004 lOkamoto et all 12005 



stream accretion (|Keres et al.ll2005l . 120091 ; IPekel et alj2009h . 

the model has proven successful for understanding global 
properties of galaxy assembly. 
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I Governato et all 120071; ICrofl et al.ll2009l; IScannapieco et al.l 
hoogt iPiontek fc Steinmet j l2009bl ; lAgertz et all l2009l h 

To date, no attempt has yielded a realistic candidate 
The dominant reason for this is the so called "angular 
momentum problem" which leads to smal l, centrally con- 
centrated discs dominated by large bulges (jNavarro fc Benj 
Il99ll ; iNavarro fc White! 1 19941 ) . Merging substructures lose 
angular momentum to the outer halo via dynamical friction, 
forcing the associated baryons to end up in the central 
parts of the proto-galaxy as a spheroid rather than a disc. 
This poses a problem for the theoretical understanding of 
extended late- type galaxies. This might in part stem from 
numerical issues: the c ommonly used Smoothed Particle 
Hydr odynamics (SPH) l|Gingold fc Monaghar] Il977l ; iLucvl 
Il977l ) technique is known to incorrectly treat boundaries, 
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hence poorly treating multiphase fluids (e.g. lAgertz et al.l 
120071 ; iRead et alj|2010l ). This can lead to artificial angular 
momentum transfer at the inter face between cold disc and 
a hot halo (|Okamoto et al 1 120051 ). 

Many proposed solutions exists to the angular momen- 
tum problem, all amounting to the same process: keep the 
gas from cooling and forming stars too efficiently in the 
merging dark matter satellites at high redshift. One nat- 
ural source is the cosmological UV background, being re- 
sponsible for reionization at z > 6 which heats the gas, pre- 
venting it to cool efficient l y into star-forming dwarf galaxies 
llThoul fc Weinberdll996l : lOuinn et all 11999 : iGnedinllioOol : 
iHoeft et al.ll2006T l. However, the impact on objects larger 
than v c i TC ~ lOkms -1 is uncle ar due to e.g. self-sh ielding 
and efficient collisional cooling l)Diikstra et al.ll2004l ). 

Gas in low-mass haloes can also be blown out by su - 
pernova driven winds (|Dekel fc Silklll986l : lEfetath iou 2000), 
hence lowering the resulting star formation efficiency (SFE), 
enriching the intergalactic medium (IGM) in the process. 
iMac Low fc Ferraral (|l999h demonstrated that while dwarf 
galaxies of mass 10 6 — 10 9 Mq efficiently can expel metals in 
supernovae-driven winds, virtuall y no mass is lost for sys - 
tems of mass > 10 7 M o (see also bubois fc TevssiedboOSh . 
The ineffi ciency in driving wind s from dwarfs was also re- 
ported bv lMarcolini et all (|2006h who attributed this to the 
extended dark matter halo and efficient metal cooling. In 
this scenario, mass-loss and IGM enrichme nt will occur due 
to ti dal and ram-pressure stripping (e.g. iMori fc Burkertl 
2000). Phenomenological models of e.g. momentum driven 
winds have proven successful in reproducing the high-z IGM 
l|Oppenheimer & Dave 2006) but it is uncertain how it reg- 
ulates star formation and in what manner the expelle d gas 
is re-accreted at later times (|Qppenheimer et al.ll201oT ). 

Various recipes of supernovae feedb ack have been de- 
veloped for numerical simulations (e.g. iNavarro fc White! 
1 19931 : iKav et all 120021 : IScannapieco et al l l2006h . and the 
methods have proven successful in removing low an- 
gular m omentum material from central parts of galax- 
ies (e.g. ISommer-Larsen et al.l 120031 : lOkamoto et alJ 120051 : 
iGovernato et al.ll2007t l. yielding more extended galaxies in 
comparison to models without feedback. However, it is un- 
clear to what extent this way of reducing star formation 
can account for disc-dominated sp i ral gal axies like the Milky 
Way. Recently IScannapieco et al.l |2009l ) demonstrated nu- 
merically, in a fully cosmological setting, how a set of 8 
Milky Way sized haloes failed to form significant discs. 
While half of the sample were early type galaxies result- 
ing from late time mergers, the other half of the sample 
had less than 20 per cent of their stellar mass in discs. 
This can be a result of the inability of the adopted feed- 
back to remove or redistribute low angular momentum ma- 
terial, but is also a strong indication that something else 
might regulate star f ormat ion at high redshift. On the same 
topic, ISawala et all (|2010h argues that modern simulations 
of dwarf galaxy format i on l|Valcke et al. 2008; S tinson et al.l 
2009; IGovernato et all [2010) all yield much larger stellar 
masses than expected from observations as well as gas-to- 
star c onversion efficiencies almost an o rder of magnitude too 
large. iDutton fc van den Boschl (|2009h found that, for SNe 
feedback to yield realistic galaxies, it must be very efficient, 
converting 25 per cent of the SN energy into outflows. If 
too strong feedback is employed, the discs can be destroyed 



by internal processes as too much material is ejected into 
the halo, preventing efficient disc reformation from cold gas, 
and possibly violat ing the upper bo unds of halo gas found in 
X-ray surveys [see iBregmanl (|2007l ) and references within] . 
In light of these studies, it is unclear if supernovae feed- 
back is the sole agent in regulating star formation. Note 
that SNe explosions can regulate star formation in galaxies 
without expelling gas, bein g a driver of galactic turbulence 
(|Mac Low fc Klessenir2004r ). 

Fundamentally, star formation is regulated by the 
availability of H2. The observed Kenn icutt-Schmidt (from 
now on K-S) relation l|Kennicuttlll998h . that relates Esfr 
to E gas , varies strongly among individual sp iral galaxies 
and can not be fit with a single power law (jBigiel et al.l 
2008). Ssfr behaves very differently for E gas greater or 
smaller than ~ 9Mq pc~ 2 , marking the tra nsition from 
atom ic to fully molecular star-forming gas (|Lerov et alJ 
2008), and is dependent on gas metallicity, dust con- 
tent, turbulence, s mall scale dumpiness an d local dis- 
sociating UV field (|McKee fc Ostrikej l2007h . The inclu- 
sion of these processes and its impact on global star 
formation in discs has recently been studied both nu- 
merically llRobertson fc KravtsovlkoOSl : iGnedin et ai1l2009l ; 
Pelupessv fc Papadopoulosll2009l ~as well as analytically (e.g. 



rapadopc 

Krumholz et al.ll2009l ) . A natural outcome of this treatment 
is an order of magnitude lower amplitude of the K-S rela - 
tion at high redshifts (z ~ 3) i|Gnedin fc KravtsovlfeoiObl ). 
This agrees well with the observ ation of damped Lya sys- 
tems (DL A: IWolfe fc Chenll2006l ) as well as Lyman Break 
Galaxies l|Rafelski et al.l l2009). This indicates that star for- 
mation can be made inefficient at high redshift, leaving gas 
for late-time star formation in a disc like environment, but 
not necessa rily by expelling gas in supernova-driven winds. 
In addition. iMurrav et alJ |2010l) argues that the disruption 
time-scale of giant molecular clouds (GMCs) due to jets, HII 
gas pressure, and radiation pressure also serves to regulate 
the SFE in galaxies. The disruption occurs well before the 
most massive stars exit the main sequence, meaning that 
supernovae in principle have little effect on GMC lifetimes. 

In this paper we investigate to what extent supernovae 
feedback and the underlying small scale star-forming physics 
can affect the formation and evolution of realistic spiral 
galaxies in a fully cosmological setting. The former effect is 
studied via well tested numerical implementations of SNII, 
SNIa feedback coupled to metal enrichment, as well as stel- 
lar mass- loss. The latter influence is achieved by considering 
different normalizations of the Schmidt-law star formation 
efficiency. We conduct a comprehensive analysis of the re- 
sulting z = discs and compare them to observational rela- 
tions. 

The paper is organized as follows. In Section[2] we de- 
scribe the numerical method used in this work, including 
the adopted feedback and star formation prescriptions. In 
Section[3l we present the cosmological initial conditions and 
discuss the free parameters of this work. Section|4] outlines 
the disc analysis and summarizes the final properties of the 
simulation suite. In Section[5] and Section[6] we present a 
detailed analysis of the impact of small-scale SFE and su- 
pernova feedback respectively. In Section[7]we compare our 
simulations to modern observations. Finally, Section[8] sum- 
marizes and discusses our conclusions. 
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2 NUMERICAL FRAMEWORK 

We use the Adapti ve Mesh Refinement (AMR) code RAM- 
SES (|Tevssierl F2002) to simulate the formation of a massive 
disc galaxy in a cosmological context including dark mat- 
ter, gas and stars. The gas dynamics is calculated using a 
second-order unsplit Godunov method, while collisionless 
particles (including stars) are evolved using the particle- 
mesh technique. The equation of state of the gas is that of 
a perfect mono-atomic gas with an adiabatic index 7 = 5/3. 
Self-gravity of the gas is calculated b y solving the Poisson 
equation using the multi-grid method l|Brandtl 19771 ) on the 
coarse grid and by the conjugate gradient method on finer 
ones. The modelling includes r ealistic recipes for star forma- 
tion (R ascra & Tcyssicr 2006), supernova feedback and en- 
richment l|Dubois fc Tevs sicr 2008). Details on these imple- 
mentations are given below. Metals are advected as a passive 
scalar and are incorporated self-consistently in the cooling 
a nd heating routine. The code adopts the cooling function 
of ISutherland fc Dopital jl993h for cooling at temperatures 
10 4 — 10 8,5 K. We extend coolin g down to 300 K using rates 
form iRosen fc Bregmanl (|l995T ). Gas metallicity is also ac- 
counted for in the cooling routin es. A UV background is c on- 
sidered using the prescription of lHaardt fc Madaul (1 19961 ). In 
order to model a subgrid gaseous equation of state, hence 
avoiding artificial gas fragmentation, the gas is given a poly- 
tropic equation of state 



T = T 



(1) 



for densities large than po. Throughout this paper we adopt 
To = 1000 K and 70 = 2.0. In this work, the polytrope den- 
sity is_set_ec[ual to the star formation threshold no. Follow- 
ing [Xgertz^etlil] l|2009l ). we adopt an initial metallicity of 
Z = 10 -4 Zq in the high-resolution region. This also serves 
as a flag for allowed regions of refinement. The refinement 
strategy is based on a quasi-Lagrangian approach, so that 
the number of particles per cell remains roughly co nstant, 
avoiding discreteness effects fe.g. iRomeo et al.1l2008h . 



2.1 Star formation 

To model th e conversion of gas into stars we adopt a 
Schmidt-law (|Schmidtll 19591 ) of the form 



Pg 



-£«— for p > po, 

iff 



(2) 



where p g is the gas density, tg — ^/37r/32Gp is the local 
free-fall time, eg is the star formation efficiency per free- 
fall time and po is the threshold for star formation. As soon 
as a cell is eligible for star formation, particles are spawned 
using a Poisson process where the stellar mass, m* , is chosen 
to be a multiple of poAa; 3 . Each formed star particles is 
treated as one stellar population with an associated initial 
mass function (IMF). This is a relevant approximation as the 
star particle masses are orders of magnitudes larger than the 
average stellar masses. We also ensure that no more than 90 
per cent of the gas in a cell is depleted by the star formation 
process. 

The po and eg parameters in Eq.[2]are, in addition to be- 
ing unconstrained physical parameters, resolution and hence 
scale dependent. There are in a sense two regimes of star for- 
mation in global simulations of disc galaxies as follows. 



(i) The ISM is resolved: At parsec-scale resolution, 
star formation occurs in their natural sites i.e. massive 
clouds such a s GMCs. Modern estima tes of star formation 
efficiencies bv lKrumholz fc Tanl (|2007l ) point towards values 
of eg = 1 — 2 per cent at densities of n ~ 10 2 — 10 5 cm" 3 . 
To only allow for star formation to take place in the 
actual physical star formation sites, hence tracing the 
formation of H2, allows for more accurate predictions of 
e.g. the Kennicutt-Schmidt star formation relation with 
less of a requirement to tune numerical parameters (see 
e.g. iGnedin et al.|[2009l ). This treatment leads to p g — > ph 2 
in Eq.[5] which is equivalent to eg being dependent on 
the environment, due in part to the local H2 fraction. 
On galactic scales, this means that the scaleheight of all 
ISM com ponents are r esolved using at least 10 resolution 
elements |Romeolll994l ). If this is not satisfied, the true disc 
stability will not be modelled accurately. This treatment is 
the goal of most simulations, but is due to the computa- 
tional load beyond the capabilities of modern simulations 
attempting to study the assembly and evolution of large 
spiral galaxies to z = 0. Isolated simulations of large spiral 
galaxies in a non-cosm ological setting have successfully 
reach ed this resolution (|Agertz et al.l 120091 ; iTasker fc Tanl 
2009), albeit with simplified physics. As the star formation 
sites become resolved, new physics becomes important e.g. 
radiative feedback i n order to accuratel y treat the lifetimes 
of GMC structures (|Murrav et al.ll2010l ). 

(ii) The ISM is under-resolved: To radially resolve a 
Milky Way like galactic disc, i.e. sampling the scale radius 
with at least 10 resolution elements, a force and hydro res- 
olution of a few 100 pc is necessary. At this resolution the 
scale height is captured with more or less one resolution el- 
ement. The true disc stability can be affected as both the 
density and velocity structure (gas and stellar dispersion) 
are influenced numerically. This still allows for the disc to 
have the correct global properties such as gas and stellar 
mass compositions, thin and thick disc, and even to develop 
realistic spiral structure. In this case a statistical star forma- 
tion recipe based on the local gas density and free-fall time 
is well motivated both theoretically and observationally. 



As we will describe in Section [37Tl we are targeting the 
latter regime of subgrid star formation and will investigate 
some of the numerical caveats related to it. At resolutions 
of several 100 pc, the eff parameter absorbs the small scale 
physics regulating star formation, allowing for a qualitative 
influence on galaxy formation. We note that there are alter- 
native formulat ion of gas to star conversion laws in the liter- 
ature (see e.g. iLerov et al.l (|2008l ) for a comprehensive sum- 
mary). However, as they are all designed to fit an observed 
relation, and all include a normalization constant similar to 
eg, we believe Eq.[2]is a representative choice for this study. 



2.2 Supernovae and stellar feedback 

The standard recipe for supernova feedback in RAMSES in- 
volves only Type II supernovae events (SNII). We have also 
implemented additional treatment of Type la events (SNIa) 
as well as mass-loss via stellar winds. Including all of these 
effects is important as a single stellar population can return 
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up to 30 — 40 per cent of its mass to the ISM during its 
lifetime. The implementations are as follows. 

2.2.1 Type II 

Type II SN events are relevant for stellar masses of 8—40 Mq 
which represents ~ 10 per cent of the mass of a stellar popu- 
lation, regardless of IMF. We assume that 10 Myr after a star 
particle is formed, 10 per cent of the star particle's mass is in- 
jected into the nearest gas cell together with a total energy of 
Esim = 10 51 (m e jecta/10 M©)erg in thermal energy. At low 
resolution, this energy would quickly radiate away jKatzl 
1 19921 ) in the dense gas, without allowi ng for an adiabatic 
expa nsion of the supernova blast-wave (|McKee fc Ostrikerl 
1977). To remedy this, we turn off cooling in cells contain- 
ing young stars to allow for the blast-wave to grow and be 
resolved by few cel ls, hence conver ting thermal energy into 
PdV work (see e.g. iGerrit scn 1997). In detail, for every star 
formation event, the inverse of the birth time, 1/tbt, is stored 
in the computational grid, overwriting any previous value. 
This passive scalar field is advected with the hydro flow (the 
conserved quantity is p gas /ibt). Gas cooling at a simulation 
time tsim is only allowed if t s lm — tbt > At g, where At Q ff is 
the cooling shut-off time-scale. Calculations of relevant time- 
scales and numerical test s usin g the SPH formalism was car- 
ried out bv lStinson et al.l (|2006f ) . Relevant time-scales are on 
the order of tens of Myr and we adopt At s — 50 Myr. 

2.2.2 Type la 

We treat S NIa and stella r mass -loss using the prescription 
outlined in iRaiteri et all (Il996l). The assumed IMF is the 
parametrization by Kroupa et al.l (|l993T l. which for a star 
particle of mass m, reads 

( nO.9,^— 1.3 



$(M) = m,A ■ 



M~ 
M~ 



2:2 
2.7 



for 0.08 < M < 0.5 M Q 
for 0.5 < M < 1 M 
for M > 1 M Q , 



where M is here the stellar mass in units of Mq and the 
normalization constant A ~ 0.3029. The adopted lower and 
upper limits are 0.08 Mq and 100 Mq respectively. At each 
simulation time-step, we calculate the mass fraction of each 
star particle ending its H and He burning phase, i.e. leaving 
the main sequence, using the fit 

log*. = a a (Z) + ai (Z) log M + a 2 (Z) (log M)\ (3) 

where t* is the lifetime of the star and Z the metallicity. 
The adopted coe fficients and r e ferenc es to the original data 
can be found in IRaiteri et all (|l996T l. Progenitors of SNIa 
are carbon plus oxygen white dwarfs that accrete mass from 
binary companions. Stellar evolution theory predict the bi- 
nary masses to be in the range of ~ 3 — 16 Mq. The number 
of SNIa events within a star particle, at a given simulation 
time with an associated timestep At, is 



$(M 2 )dM 2 , 



(4) 



where m t and m t +At is the mass interval of stars ending 
their life during the computational timestep. ^(M 2 ) is the 
IMF of the secondary star, i.e. 



4>(M 2 ) = m,A' 



A/sup 



M inf 



M 2 
Mb 



where Mb is the mass of the binary, Minf = max(2M2, 3 Mq) 
and M Bup = M 2 + 8M Q . The constant A' = 0.16 A, 
and is a calibrated value for SN Ia events in our Galaxy 
(van den Bergh fc McClurd 1 19941 ). Each explosion is as- 
sumed to release 10 51 erg (released as thermal energy in 
the nearest gas cell) and 0.76 Mq of meta l enriched material 
(0.13 M s of ie O and 0.63 M Q of 56 Fe) (|Thielemann et all 
Il986l ). 



2.2.3 Stellar mass-loss 

For each time-step At, and star particle, we calculate the 
average stellar mass, (M), exiting the main sequence using 
Eq.[3] The mass-loss during the At time-span is calculated 
using the best fit initial- final mass relation of iKalirai et al.l 
i|2008l ) : 

M wind = 0.891 - 0.394/(M) . (6) 
At each time interval At, the total mass-loss in winds is 

A^tot.wind = fmrn^M^nd, (7) 

where 



fn 



fm-t + At 
J mt 



$(M)dM. 



(8) 



The lost stellar mass enters the gaseous mass in the nearest 
cell and the gas metallicity is updated consistently with the 
star particle's metallicity. 



3 INITIAL CONDITIONS AND SIMULATION 
SUITE 

The initial conditions used in this work are a subset the 
Silver River simulation suite (Potter et al. in preparation), 
aimed to study the pure dark matter assembly history of 
a Milky Ways size h alo in much greater d etail than here. 
We adopt a WMAP5 l|Komatsu et al.ll2009l ) compatible cos- 
mology, i.e. a ACDM Universe with Q A = 0.73, Q m = 0.27, 
fi b = 0.045, CT 8 = 0.8 and H = 70kms _1 Mpc" 1 . 

The upcoming work by Potter et al. will present the 
details concerning the initial condition generation. Briefly, 
a pure dark matter simulation was performed using a sim- 
ulation cube of size Lb ox = 179 Mpc. At z — 0, a halo of 
mass M 20 o,c w 9.7 x 10 11 M was selected for re-simulation 
at high resolution, and traced back to the initial conditions 
at z = 133. M2oo,c is the virial mass of the halo, defined 
as the mass enclosed in a sphere with mean density 200 
times the critical value. The corresponding virial radius is 
?~200,c = 205 kpc. By using a definition based on 200 times 
the background density we obtain M200.bg = 1-25 x 10 12 Mq 
and r2oo,b g = 340 kpc. When baryons are included in the 
simulations, the final total halo mass remains roughly the 
same. 

The halo has a quiet merger history, i.e. it undergoes 
no major merger after z = 1, which favours the formation of 
a late-type galaxy. A nested hierarchy of initial conditions 
for the dark matter and baryons was generated using the 
GRAFIC+40 code, where we allow for the high resolution 



(5) 



1 http://grafic.sourceforge.net/ 
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particles to extended to 3 virial radii from the centre of the 
halo at z = 0. This avoids mixing of different mass dark mat- 
ter particles in the inner parts of the domain. In this work, 
we focus on two sets of resolutions from the Silver River 
suite, referred to as SR5 and SR6. The simulations are iden- 
tical apart form the number of dark matter particles, and 
hence particle mass, as well as maximal AMR refinement. In 
SR6 the dark matter particle mass is moM = 2.5 x 10 6 Mq 
and in SR5 moM = 3.2 x 10 5 M©. The mesh is refined if 
a cell contains more than eight dark matter particles, and 
similar criterion is employed for the baryonic component. 
At the maximum level of refinement, the simulations reach 
a physical resolution of Ax = 170 pc and Ax = 340 pc in 
SR5 and SR6, respectively. 

3.1 The free parameters 

The goal of this work is to study how the characteristics of 
disc galaxies change when standard numerical parameters 
governing star formation are modified. Following the discus- 
sion in Section[TJ we consider star formation regulation in 
two different ways: small scale (~ 100 pc) physics such as 
H2 abundance, GMC turbulence, metallicity, radiative ef- 
fects etc., or via energy injection from supernovae explosions 
leading to gas expulsion in galactic winds. The first mech- 
anism is modelled by varying the Schmidt-law (Eq.[2|) SFE, 
eff, which acts on a cell- by-cell basis. The latter is studied 
by increasing the injected SNII energy, JSsnii- In addition 
we study the impact of the star formation threshold no, but 
in lesser detail. 

The traditional way of treating star formation in sim- 
ulations of galaxy formatio n (e.g. iGovernato et all 120071 ; 
iPiontek fc Steinmet j|2009bh is to tune the SFE parameter 
using an isolated disc model to match the observe d K-S rela- 
tion, m ost commonly the fitting formula given bv lKennicuttl 
( 1998) of z = galaxies. In addition, the recipe for energy 
injection via supernovae and its efficiency is tuned simulta- 
neously. These parameters are then used in fully cosmolog- 
ical simulation of galaxy formation. This scheme assumes 
that supernova explosions are the main sources of star for- 
mation regulation at high redshift. As argued in Section]]] 
the numerically assumed constant efficiency is strongly red- 
shift dependent and a z = tuning is likely to over-predict 
star formation in more metal poor environment at higher 
redshift. We treat eg as a free, but constant, parameter 
and adopt eg = 1,2 or 5 per cent in the fully cosmologi- 
cal context. These values are i n agre ement with GMC es- 
timates from iKrumholz fc Tan! ((2007). As we will demon- 
strate below, lower values than what traditionally is adopted 
is preferred in order to form late-type galaxies. Note that 
eg ~ 2c* , where c* is t he e fficiency parameter used i n e.g. 
IGovernato et al.1 (120071 ) and IScannapieco et al.l (|2009h (de- 
fined via idyn = l/s/4irGp). Values of c» = 0.05 — 0.1 are 
commonly employed i.e. a few times, up to an order of mag- 
nitude larger than what we consider here. 

The standard SNII feedback described in Section [2.2.1l 
is the baseline feedback in all of our simulations. In a subset 
of simulations we add the additional stellar mass-loss and 
SNIa treatment. The high-efficiency simulation (eg = 5 per 
cent) is used as a template for the impact on the injected 
feedback SNII energy, which we set to -Esnii = 1,2 and 
5 x 10 51 erg. Using energies that are several times larger than 



the canonical 10 erg might be perceived as unrealistic, but 
we believe it is illustrative to study the extreme cases of this 
type of feedback. In addition, the amount of SNII energy 
dissipated in cooling, after the shut-off time has passed, is 
complicated to measure. As a control set we also run the 
simulations without feedback, both with and without metal 
enrichment. 

The philosophy of the star formation threshold is as 
follows. In reality stars form in molecular clouds of average 
densities of n > 10 2 cm -3 . Imposing a threshold of this mag- 
nitude would require a resolution on the order of parsecs to 
resolve the formation of the star-forming clouds, something 
that is beyond the scope of fully cosmological hydr o+A^- 
body simulations today (but see iGnedin et all l2009h . We 
adopt no = 0.1 and 1cm" 3 for each setting of eg, but the 
appropriate choice is fundamentally tied to resolution and 
can lead to spurious results. The ISM has been shown to 
be represented by a logno rmal density probability distribu- 
tion function (PDF) (e.g. iKravtsovl 12003 ; IWada fc Normanl 
l2007h . or even a super position of several lognormal ly dis- 
tributed ISM phases (|Robertson fc Kravtsovl 120081 ). The 
amount of gas eligible for star formation is represented by 
the high density part of the PDF which in turn is a function 
of total disc gas mass and turbulence. A density threshold 
should be picked to allow for the high-density star-forming 
part of the PDF to be well resolved or at least contains, 
given an adopted numerical resolution, a converged amount 
of star-forming mass. If not, then chosen threshold will af- 
fect the numerical efficiency for global star formation. We 
will demonstrate this effect below. 

In summary, the varied constants of interest is here the 
star formation threshold (no), the star formation efficiency 
per free-fall time (eg) and the form of supernova feedback 
and injected energy (J5snii). Our main focus is the impact 
of these parameters at the SR6 level of resolution, and we 
present a brief resolution study in AppendixlAl We summa- 
rize our complete test suite in Table [1] 



4 THE DISCS 

In this work, we focus primarily on the disc properties in 
the SR6 simulations at z = 0. Details of the satellite galax- 
ies and halo properties will be considered in a future work. 
FigsQ] and [2] show projected face-on and edge-on stellar and 
gas density maps at z — for the galaxies in the star forma- 
tion efficiency and feedback test suite, respectively. The discs 
show a wide range of spiral galaxy morphologies, and will 
return to this point in Section [5~4l We decompose the result- 
ing stellar discs into a bulge, bar and disc component and 
fit these simultaneously to the stellar surface density pro- 
file. The latter is calculated using all stars out to a height of 
\z\ = 2.5 kpc. For the bulge and disc component we assume 
exponential profiles, i.e. 



E(r) = £ exp(-r/r d ), 



(9) 



where we fit for Eo and the scale radius r&. The bar compo- 
nent is modelled using a simple Gaussian, 



E bar (r) 



A 



■ exp 



(r - r f 
2a 2 



(10) 



where we fit for the width a, the central point ro and am- 
plitude Aq. We consider this a conservative estimate of the 
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Table 1. Summary of the numerical parameters. The simulations use a maximum physical cell resolution of Ax = 340 pc (SR6) or 
Ax = 170 pc (SR5), and the high-resolution region is occupied with dark matter particles of mass mrjM = 2.5 X 10 6 Mg (SR6) or 
moil = 3.25 X 10 s Mg (SR5). All simulations use delayed cooling in regions of young stars, unless specified. When SNII feedback is used 
-E'SNII = 10 51 erg, unless other values are indicated. 



Run 




Feedback 


Star formation threshold, no (cm 3 ) 


SR6-n01el 


1 


SNII 


0.1 


SR6-n01e2 


2 


SNII 


0.1 


SR6-n01e5 


5 


SNII 


0.1 


brlD-nUlellVlij 


1 


oJNll, olNla, mass- loss 


0.1 


SR6-n01e2ML 


2 


SNII, SNIa, mass-loss 


0.1 


O f* J~\ 1 f TV J T 

SR6-n01c5ML 


5 


SNII, SNIa, mass-loss 


0.1 


SR6-nlel 


1 


SNII 


1 


SR6-nlc2 


2 


SNII 


1 


SR6-nlc5 


5 


SNII 


1 


SR6-nlelML 


1 


SNII, SNIa, mass-loss 


1 


SR6-nle2ML 


2 


SNII, SNIa, mass-loss 


1 


SR6-n01elNFB 


1 


No feedback, Esnii = 


0.1 


SR6-n01e5NFB 


5 


No feedback, Esnii = 


0.1 


SR6-n01e5NFBmet 


1 


No feedback but metal enrichment 


0.1 


SR6-n01e5SN2 


5 


SNII, BgNn = 2 X 10 51 erg 


0.1 


SR6-n01e5SN5 


5 


SNII, E sml = 5 x 10 51 erg 


0.1 


SR5-nlclML 


1 


SNII, SNIa, mass-loss 


1 



Table 2. Summary of disc characteristics at z = 0. The mass of the components are obtained by fitting the stellar surface density (see 
text), and are in units of 10 10 M . Note that we consider all gas phases for the gas mass and all stars for the stellar masses. (1) Fitted 
scalclcngth of stellar disc. Large uncertainties exist for > 10 kpc as the stellar discs are small and feature almost flat stellar surface 
density profiles. (2) / gas = Mdi sc , g /(Mdisc,tot + Mb u i g o). (3) Total measured specific angular momentum of the baryons in the disc and 
bulge in units of kms -1 kpc. 



Run 




Afdiac.g 


-Mbulge.s 


r d (kpc) (1) 


/gas (2) 


B/D 


B/T 


Jbar (3) 


SR6-n01el 


8.6 


1.6 


2.0 


3.8 


0.13 


0.23 


0.19 


1920 


SR6-n01e2 


7.4 


1.3 


4.6 


7.6 


0.10 


0.62 


0.38 


1655 


SR6-n01e5 


5.6 


0.72 


7.0 


~ 15.0 


0.05 


1.25 


0.56 


1305 


SR6-n01clML 


8.0 


2.3 


2.2 


5.0 


0.18 


0.27 


0.21 


1960 


SR6-n01c2ML 


8.1 


1.6 


3.8 


5.0 


0.12 


0.47 


0.32 


1718 


SR6-n01e5ML 


5.5 


0.93 


7.2 


~ 15.0 


0.07 


1.30 


0.57 


1464 


SR6-nlcl 


6.6 


3.3 


2.9 


2.7 


0.26 


0.44 


0.31 


1594 


SR6-nlc2 


6.4 


2.4 


4.3 


2.5 


0.18 


0.67 


0.40 


1804 


SR6-nlc5 


6.0 


2.1 


5.2 


2.7 


0.16 


0.87 


0.46 


1643 


SR6-nlclML 


6.5 


3.6 


2.7 


2.7 


0.28 


0.42 


0.29 


1618 


SR6-nle2ML 


6.3 


2.9 


4.3 


2.7 


0.21 


0.68 


0.41 


1281 


SR6-n01elNFB 


7.8 


1.6 


3.1 


4.0 


0.13 


0.40 


0.28 


1938 


SR6-n01e5NFB 


5.8 


0.62 


6.5 


~ 15.0 


0.05 


1.12 


0.53 


1394 


SR6-n01e5NFBmct 


5.6 


0.56 


6.6 


10.0 


0.05 


1.18 


0.54 


1430 


SR6-n01e5SN2 


5.6 


0.55 


6.5 


4.5 


0.05 


1.16 


0.54 


1266 


SR6-n01e5SN5 


6.6 


0.81 


2.3 


2.8 


0.09 


0.35 


0.26 


1387 


SR5-nlelML 


9.0 


2.0 


2.2 


2.8 


0.18 


0.24 


0.20 
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Figure 1. Projected face-on and edge-on surface density maps of the stars (top) and gas (bottom) of the 2 = discs, where each panel is 
60 kpc across. As the SFE is lowered and mass-loss is employed, spiral structure becomes more pronounced due to a less massive bulge. 
The Hubble type of the disc changes from an early type (SO or Sa) disc, to a late- type spiral galaxy (Sb or Sbc) as we decrease eg . 




Figure 2. Projected face-on and edge-on surface density maps of the stars (top) and gas (bottom) of the z = discs where each panel is 
60 kpc across. From left to right, the feedback strength is EsNII = 0, 1, 2 and 5 X 10 51 erg. As the injected amount of energy is increased, 
the gas component becomes more distorted. 



bar mass as a Gaussian contribution falls off towards the 
centre of the disc, leaving more mass to be accounted for 
by the bulge. An example of the fitting procedure can be 
seen in Fig. [3] The necessity of a separate bar component is 
here clearly illustrated. In the more bulge-dominated cases 



the bar amplitude is decreased considerably, owing to the 
weaker disc self-gravity. 

The bulge mass, A/buige — 27rSbuige?*buige; is obtained by 
integrating Eq.[5] A similar relation holds for the disc, where 
we also include the stellar disc mass past the break radius 
in the quoted disc stellar mass, Mdisc.s, but we only use the 
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Figure 3. Example of a multi-component fit (red dashed line) 
to the stellar surface density (black solid line) of the z = disc 
in nOlelML. The fit is composed of a bulge (red dotted line), a 
bar (blue dotted line) and a disc (black dotted line) component. 
The green dashed line shows the bulge and bar surface density 
contribution. 



data within the break when fitting (see Fig-[3jl . The bar mass 
is simply the integrated mass from Eg. 1101 and we consider 
the bar as a part of the disc component and include it in the 
quoted Mdi sc ,s- Doing this or not modifies the disc mass only 
slightly, especially in bulge-dominated galaxies. The gas is 
treated as a single component, and we simply consider the 
mass within r = 15 kpc and \z\ = 2.5 kpc as the gaseous disc 
mass, Mdisc, g - We consider only the stars when calculating 
the bulge-to-disc (B/D) and bulge-to-total (B/T) ratios. All 
measured and derived quantities are summarized in Table [2] 
We note that this method of defining galactic components 
in si mulations, as well as others e.g. via angular m omen- 



turn (Oka moto et al.ll2005l ; IScannapieco et al. |2009| ). carry 
uncertainties. 

Our simulated discs span a large range of characteris- 
tics: stellar disc masses are in the range Mdiscs = 5 — 9 x 
10 10 M , bulge masses of M bu i gc , s = 2 - 7 x 10 10 M©, B/D 
~ 0.23—1.2 and gas fractions / g = 0.05—0.28. The scale radii 
of the discs, rd, vary from typically 4 — 5 kpc to > 10 kpc 
in the bulge-dominated systems. As we demonstrate below, 
extended disc galaxies of Sb, or even Sbc type, form only 
when star formation is numerically resolved in the whole 
disc (no = 0.1 cm -3 ), and a low efficiency of eg ~ 1 per cent 
(or very strong feedback) is adopted. At larger efficiencies 
we observe how the B/D ratios increase, the discs are less 
extended, the rotational velocities peak at very large values 
and the spiral patterns become more tightly wound and less 
pronounced. This indicates a shift towards early type discs 
like Sa or even SO. 



5 EFFECT OF STAR FORMATION 
PARAMETERS 

In this section we study the influence of star formation pa- 
rameters, i.e. in essence the small scale physics, on disc prop- 
erties at z = 0. The resulting stellar surface densities (£ a ), 
gas surface densities (E gas ) and rotational velocities mea- 



sured from the gas (w ro t) are presented in Figs [4] and [5] for 
the first 11 simulations in the SR6 suite at z = (see Ta- 
ble. [TJ. 



5.1 The star formation density threshold, no 

We start by focusing on the data presented in Fig. [4] By 
keeping eg fixed to 1 per cent, while varying no, we ob- 
serve a strong change in the ability to form stars at large 
radii. The galaxies adopting a large threshold have a more 
concentrated distribution of stars, smaller stellar disc scale- 
lengths as well as larger £ gas at all radii. The scalelengths 
are r& > 4, kpc for no = 0.1cm -3 , but only rd ~ 2.5 kpc 
for no = 1 cm -3 . The latter values are on the low side when 
compa red to observations of late-type spiral s at this mass 
range (|Courteau et alj 1 19961 : ICourteaul Il997l ; iGnedin et al.l 
120071 ) . The systematically lower E s signals an under-resolved 
or missing' star formation throughout the disc: the average 
physical gas density does not efficiently cross the targeted 
no, even at intermediate radii. This is also reflected in the 
gas fractions of ~ 25 per cent, which is much lar ger than 
observed average value s for galaxies of this size (|Garnettl 
120021 : Izhang et al]|2009l ). The rotational velocities are large 
at smaller radii for no = 1cm -3 , regardless of choice of eg. 
Naively one would expect this numerically induced star for- 
mation deficiency to alleviate the angular momentum loss 
at high redshift, hence forming a less concentrated galaxy. 
However, due to secular evolution in the disc, this is not the 
case: disc instabilities drive gaseous flows to the centre of the 
disc where, as the gas crosses the correct threshold, star for- 
mation can proceed. We conclude, that given our numerical 
resolution and simulated system, no =0.1 cm -3 yields more 
realistic (average) disc galaxies when compared to observa- 
tions. We discuss this numerical effect and its relationship 
to the adopted mesh resolution further in AppendixlBl Note 
that this is not a fundamental result of galaxy formation, but 
serves only as tuning given our numerical resolution and is 
a basis for the subsequent tests. A discussed in Section l2~Tl 
no should be increased as the resolution is increased. 



5.2 The star formation efficiency, eg 

We now turn to the data presented in Fig.[5] where we keep 
the threshold fixed at no = 0.1cm -3 and adopt eg = 1, 2 or 
5 per cent. As eg is increased, S sta r increases at small radii 
i.e. the bulge mass increases, singalling a lower disc angular 
momentum. The bulge to disc ratio increases from B/D = 
0.25 to 1.25 as es increases from 1 to 5 per cent. £ gas roughly 
follows a 1/r- profile and the magnitude is lowered at all 
radii by approximately the relative change in efficiency. The 
stellar disc is less extended and the exponential scale-length 
increases for larger efficiencies (see Table. [2}. For eg — 1 per 
cent, the disc scalelength is measured to be rd ~ 4 — 5 kpc, in 
good agreement wit h observed average values from the SDSS 
jGnedin et ai1l2007n , while for eg = 5 per cent, rd ~ 15 kpc 
is a > 2a outlier. Large uncertainties exist for rd > 10 kpc 
as the stellar discs are small and feature almost flat stellar 
surface density profiles. 

At large radii in all simulations, rd shifts to ~ 
2 kpc. Disk breaking i s a well observed phenomenon 
l|Pohlen fc Truiillol [2006) correlated with a dip in the 



The formation of disc galaxies 9 




Figure 4. The effect of the Schmidt-law star formation density threshold. The panels show the stellar surface density (left), gas surface 
density (middle) and rotational velocity measured from the gas (right). We consider all material within a height of z < 2.5 kpc for all 
components. The star formation efficiency is eg = 1 per cent in all simulations. The different colours are described in the first panel, and 
a dashed line indicates that we use the extended feedback model (see text). 
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Figure 5. The effect of the Schmidt-law star formation efficiency. The panels show the stellar surface density (left), gas surface density 
(middle) and rotational velocity measured from the gas (right). We consider all material within a height of \z\ < 2.5 kpc for all components. 
The Schmidt-law density threshold is n = 0.1cm -3 in all simulations. The different colours are described in the first row, and a dashed 
line indicates that we use the extended feedback model (see text). 



star formatio n rate , and has been studied numerically by 
iRoskar et alj ([2008). As larger star formation efficiencies 
lead to less extended discs, the disc breaks occurs at smaller 
radii: rb rca k ~ 16, 14, 10 kpc for eg = 1,2 and 5 per cent, 
respectively. The breaks can also be seen from the average 
stellar ages (t*), shown in Fig. [6] The central parts of the 
discs generally consist of older stars formed at z > 1, and 
(i*) decreases with larger radii, reaching {£») w 6Gyr. Past 
the disc break, older stars appear which in part can be at- 
tributed to stellar migration as well as pollution by old halo 
stars with (t») w 11 — 12Gyr. We note that (i„) flattens or 
even declines towards the centre of the discs as eg is lowered. 
This is due to secular evolution: as spiral structure is more 
pronounced (as in nOlel), gas is transported towards the 
centre more efficiently and late-time star formation occurs; 
see Fig.[n 

The efficiency has a strong impact on the rotational 
velocity. The rotation curve in the n01e5 simulation features 
a strong peak in the inner parts of the disc. As eg is lowered, 
B/D decreases and the velocity profile flattens. Only when 
eff < 2 per cent can a flat rotational velocity profile be 



produced! In nOlel, the rotational velocity reaches v TO t ~ 
275 km s -1 and stays roughly flat. For n01e2 and n01e5, w ro t 
peaks at « 310 and 360 km s -1 , but converges at 275 km s -1 
close to r = 20 kpc. 

In the left-hand panel of Fig. [7] we plot the circular 
velocities, i.e. v c (r) = yj GM{< r) Jr, for the dark matter, 
disc + bulge baryons and total mass of the galaxies. We 
note that the i> c -profiles are well traced by the cold gas ro- 
tation curve (the stellar rotational velocities are lower than 
what can be expected from v c due to a larger velocity dis- 
persion). While all simulations converge at large radii, and 
have an equal dark matter and baryon contribution within 
r ~ 17 kpc, the mass distribution (and angular momentum 
distribution) differs dramatically leading to large difference 
in circular velocities. As we will demonstrate in the next 
section, a majority of the mass within the bulge component 
originates from the intense star formation epoch at z ~ 2 — 3 
where the value of eg matters the most. We also observe a 
significantly enhanced dark matter contraction at large eg. 
This effect in our simulation suite, and its relevance for di- 
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Figure 6. Average stellar ages as a function of radius. The colours 
represent different star formation efficiencies, eg= 1 (black), 2 
(red) and 5 per cent (blue). Dashed lines mark the measured 
stellar surface density breaks. 



rect dark matter detection, has recently been analysed by 
iPato et alJ (|2010l ). 



5.3 Star formation histories 

The left-hand panel of Fig. [8] shows the star formation histo- 
ries for all stars belonging to the discs in nOlel, n01e2 and 
n01e5 at z = 0. The average SFR at the current epoch is 
~ 3 — 4 Mq yr~ in all simulation, regardless of numerical 
setting. Moreover, the star formation history during the qui- 
escent phase of disc evolution, i.e. after z ~ 1, is relatively 
flat and roughly the same in all simulations. Significant dif- 
ferences occur at intense epochs of star formation, especially 
at z — 3 where the proto-disc is assembled via cold streams, 
satellite mergers and gas accretion from the hot halo, as 
seen in the simulation snapshot in Fig.|^|. Here the SFR 
peak changes dramatically from ~ 43 MQyr -1 in n01e5 to 
~ 23 M Q yr" 1 in nOlel. 

In n01e5, stars form efficiently everywhere, even in 
satellites. The gas is quickly consumed locally during the 
high-redshift assembly, and merging systems lose angular 
momentum to the dark halo, ending up in the central part 
of the galaxy. Accretion via cold streams and from the hot 
gase ous halo will still s upply the galaxy w ith unprocessed 
gas (|Dekel et al.1 120091 ; I Agertz et ail 120091 ') . but now in a 
more bulge-dominated environment. In nOlel, star forma- 
tion is less efficient and a significant portion of the mass in 



2 The image is an RGB composite image where red is tempera- 
ture, green is metals and blue is density. Each quantity is a mass 
weighted average along the line of sight. For each image pixel, we 
calculate the RGB triplet as 

/log(T/T m ) log(Z/Z m ) log(p/p m r 



(R, G, B) = 255 



V AT 



where log{T m [K], Z m [Z@], pm I cm 
{AT, AZ, Ap} = {2, 2, 6}. 



AZ 



Ap 

{4.1,-3,-4.8} and 



the merging clumps is in gaseous form. This material is lost 
to the hot gaseous halo via ram-pressure and tidal strip- 
ping, or expelled during SNe events, which later cools down 
to join the disc. Hence, the material that is not consumed 
by star formation at z — 3 is processed at a later epoch, 
closer to z = 1 — 2 (see Fig.|5Jl, but now in a more disc- like, 
higher angular momentum configuration. In fact, the trend 
at z = 3 is reversed at this later epoch, and the largest SFR 
is found for eg = 1 per cent. These two modes of star forma- 
tion are related to the cl assical angular momentum problem 
(INavarro fc Whitdll994M . and they lead to fundamentally 
different modes of disc assembly and morphology. 

A confirmation of the above discussion is shown in 
Fig. 1101 where contours of formed stellar masses are outlined 
in the star formation time vs. disc radius plane. Note that 
this mass refers to all the stars at z = contained in the disc, 
and is hence the sum of the stars formed in merging satellites 
as well as in situ. While the formed stellar mass in nOlel is 
smoothly distributed in a roughly exponential profile across 
the disc at all times, without a clear sign of extreme star for- 
mation bursts, the n01e5 simulation shows a strong central 
concentration of stars formed at t = 11.5 Gyr (z ~ 3). This 
analysis confirms the notion of efficient star-forming satel- 
lites loosing angular momentum and being dragged into the 
central parts of the galaxy. 

5.4 Hubble types 

Fig.[T] shows mass weighted projections of the stellar and 
gas surface densities for nOlelML, nOlel, n01e2ML, n01e2, 
n01e5ML and n01e5. The gaseous discs are thin and ex- 
tended in all simulations, and are surrounded by a warped 
layer of cold/war m gas, probably associate d with misaligned 
accretion events (|Shen h, Sellwood| [2006). A hot gaseous 
halo surrounds the discs, and a temperature projection (not 
shown) reveals an extended disc-halo interface of warm/hot 
gas. We will explore this in future work. 

As discussed in Section [5~2l we find a very strong trend 
in disc and bulge mass with increasing star formation effi- 
ciency. The nOlel simulation feature a 8.6 x 10 10 Mq stellar 
disc with a 2 x 10 10 M bulge, hence B/D~ 1/4. In n01e5 
the disc is 35 per cent less massive and the bulge 3.5 times 
more massive with B/D~ 1.25. Roughly the same scaling 
holds when including additional SNIa feedback and stellar 
mass-loss. 

All discs show spiral pattern in the gas component with 
a larger amplitude in the more gas-rich discs, having lower 
eg. We also observed spiral structure in the stellar compo- 
nent, which is the most pronounced in nOlel and nOlelML 
and n01e2ML simulations. As eg is increased, B/D increases 
and the spirals arms become more tightly wound as marginal 
gravitational instabilities can no longer excite pronounced 
open spiral arm structure. All discs feature a stellar bar, 
and viewed edge-on, we observe how the inner stellar dis- 
tribution flattens as eg is decreased. The flattened central 
parts of nOlel and nOlelML, and the fact that the bulge 
is well fitted using an expon ential, is indicative of a bulge 
formed via secular processes (iKormendv &: Kennicuttll2004l ') 
e.g. bar buckling (|Debattista et al.l 20061 ). The gaseous bar 
strengthens at lower eg, and in nOlel and nOlelML gas is 
transported towards the disc centre, triggering star forma- 
tion. In these simulations, close to 50 per cent of the stars 
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Figure 7. Circular velocities, v c = \J GM(< r)/r, of the stars (solid), dark matter (dotted) and total mass (dashed). (Left) The effect of 
the star formation efficiency: eg= 1 (black), 2 (red) and 5 per cent (blue). (Right) The effect of supernovae feedback strengths (adopting 
eff= 5 per cent): Esnii = (blue), 10 51 erg (black), 2 X 10 51 erg (green) and 5 X 10 51 erg (red). 
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Figure 8. (Left) Star formation histories using different values of eg. At z < 1, all simulations regulate to similar SFRs, ending up 
at ~ 3 — 4 Mq at z = 0, regardless of star formation parameters. Using a high efficiency leads to central galaxies and dwarfs burning 
their fuel quickly at high redshift during galaxy assembly, resulting in excess angular momentum loss and a prominent central spheroid. 
A lower efficiency avoids this issue, leaving more gas left for star formation at lower redshifts in a more disc-like configuration. (Right) 
Star formation histories for a set of simulations of increasing supernovae feedback strength (-Esnii). We note that the large SFR peak at 
z = 3 is only lowered when a very large amount of energy is injected into the ISM. 



associated with the bulge formed in situ of the disc at z < 1, 
and only ~ 25 per cent formed at the intense star formation 
peak at z ~ 3. This indicates that a significant portion of 
the flattened bulge has formed via secular evolution, lead- 
ing to a pseudo-bulge. This is in stark contrast to the bulge 
formation epoch seen in the central parts of the n01e5 sim- 
ulations (right-hand panel of Fig. l 10[l . where essentially all 
bulge stars form at z ~ 3. 

IWeinzirl et al.l |2009l ) (see also lLaurikainen et"aT]|201fi ). 

recentl y analysed 182 Jf-ba nd images from the OSUBSGS 
survey (|Eskridge et al.ll2002l ) to obtain B/D and B/T values 
across the Hubble sequence. Comparing their sample aver- 



ages (see e.g. their fig. 14) to our set of simulations (B/D 
in Table[2]) suggests that the final disc in n01e5 is of SO/a 
type, in n01e2 it is of Sa/Sab type and in both nOlel and 
nOlelML it is of Sb/Sbc type. We consider this agreement 
only as indicative as ea ch Hubble type spans a wi de range of 
B/D and B/T values. iGraham fc Worlevl (|2008l ) presented 
B/D and B/T flux ratios using a sample of over 400 galax- 
ies observed in the K band. Their B/D estimates for differ- 
ent Hubble types confirm the classification of our simulated 
discs. There is no doubt that we are measuring a transfor- 
mation along the Hubble sequence. 
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Figure 9. A large scale view of the assembling spiral galaxy from the SR5 simulation at z ~ 3; the most intense epoch of star formation 
for this system. The i?G-B-image 2 shows the gas component using temperature (red), metals (green) and density (blue). We can clearly 
distinguish accretion via streams of cold pristine gas (in blue) penetrating the shock- heated gas (in red), reaching the heart of the halo. 
Dwarf galaxies outside of the large gaseous halo are surrounded by puffy distributions of enriched gas originating from stellar outflows. 
Gas is efficiently lost via tidal and ram-pressure stripping as the dwarfs interact with the main galaxy and its hot gaseous halo. The 
distance measure is in physical units. 




r [kpc] r [kpc] 

Figure 10. Formation time of stars and their radial distribution for the discs in nOlel (left) and n01c5 (right) at z = 0. The contours 
trace regions of binned mass using bins of size At = 0.5 Gyr and Ax = 0.5 kpc. The contour lines trace, from thin lines with light shades, 
to thick lines with dark shades, the formed stellar masses from log(M*) = 6.5 to log(Af») = 9.5 in steps of 0.25 dex. While the formed 
stellar mass in nOlel is smoothly distributed across the disc at all times, the n01e5 simulation shows a strong central concentration of 
stars formed at t = 11.5 Gyr (z ~ 3). 



6 EFFECT OF SUPERNOVA FEEDBACK 

As we have demonstrated in the previous sections, a high 
SFE overproduces the central stellar mass of the galaxy. 
The inclusion of additional SNIa feedback and stellar mass- 
loss did not drastically change the galaxy properties, even 
though differences can be seen in Fig.[T] (the discs in nOlel 
and n01e2 feature much stronger spiral structure) and Fig. [4] 
and more late time star formation is made possible (see 
Section | 7.2[l. The eff e ct of stellar mass-loss was studied by 
iMartig fc Bournaudl (|2010l ) who found a stronger effect on 
the bulge mass in a similar setting, perhaps due to imple- 
mentation differences. 

In Figs [7]and[lT]we present v c and E s for the n01e5 sim- 
ulation, but with different amounts of injected SNII energies; 



see Table[2] Note that we still enrich the ISM with metals 
in the simulation with E'snii = 0. Without this, the effect of 
metal cooling will not be present in all simulations. We find 
that the bulge mass is lowered as we increase Ssnii, but only 
for a very large injected value of 5x 10 51 erg can the disc rota- 
tion curve peak at a reasonable v c < 300 kms -1 , resembling 
that of the nOlel simulation. As for the standard feedback 
runs in the previous section, the dark matter halo is more 
contracted as star formation is less regulated. The difference 
in v c at r = 20 kpc between the n01e5SN5 and the other 
simulations, corresponding to a few 10 10 Mq , is due to the 
expelled gas during galaxy assembly which can be accounted 
for in the more massive gas halo. This effect can also be seen 
in E s as the central values are decreased, the discs scale ra- 
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Figure 11. Stellar surface densities of the 2 = discs in the 
supernova feedback test suite. As the SNII feedback energy in- 
put is increased, the disc becomes more extended and the bulge 
component less massive. 



dius decreases and the break radius is shifted to larger radii. 
As seen in Tabled a massive disc still forms. The effect 
on the star formation histories are shown in the right-hand 
panel of Fig. [8] we find no significant difference among the 
simulations, apart for the very energetic n01e5SN5 simula- 
tion. The SFH now resembles that of nOlel where the z = 3 
amplitude is lowered to « 2OM0yr -1 and more gas is left 
to form stars in a disc-like environment at z ~ f — 2. 

The projected gas density and stellar maps were shown 
in Fig. [2] While the standard feedback simulations shown 
in Fig.[T] showed a clear Hubble sequence of open to tightly 
wound spiral structure as B/D was lowered, this is not the 
case for the feedback test suite. In n01e5SN5, the gaseous 
disc is heavily distorted, warped and puffed up by the large 
SNII energy injections. Star formation is here very different 
compared to nOlel as stars form in filaments and shells from 
SNe explosions rather than in gas-rich spiral arms. 

The effect of metal co oling is not always accounted f or in 
cosmological simulations. iPiontek fc Steinmet j l|2009al ) in- 
cluded this effect and reported on difficulties in suppress- 
ing the initial high-z peak, even with sophisticated feedback 
models. In our simulations metal cooling is roughly counter- 
acted by the standard SNII feedback. If metal enrichment is 
turned off together with the feedback, we do not find a signif- 
icant modification to our discs. For example, the nOlelNFB 
simulation shows a surprisingly successful set of characteris- 
tics when compared to nOlel (see Table[5J). As a zero metal- 
licity gas cools inefficiently below 10 4 K, as well as in the 
range 10 5 K to 10 7 K, the nOlelNFB disc essentially behaves 
as higher metallicity counterpart but with SNII heating bal- 
ancing cooling. This is the philosophy behind su bgrid mul- 
tiphase models (e.g. ISpringel fc HernquistJ [2003) in which 
feedback is implicitly treated as a stiff gas equation of state. 
Note that our polytropic EOS is slightly stiffer than what is 
usually adopted (7 = 2 instead of 7 = 5/3). 



7 RELATIONSHIP TO OBSERVATIONS 

7.1 Angular momentum of the baryons 

For each galaxy we calculate the cumulative specific angular 
momentum vector, defined as 



Jbar(< r) 



M(< r) 



N(<r) 

S~] rmXi X Vi 



(11) 



including all bulge and disc baryons. Here Xi and Vi are 
positions and velocities of the gas cells and star particles of 
the N elements within a radius r encapsulating the mass 
M(< r). The resulting jb ar = |j' bar | are presented in Table 

El 



Using the s ample of Courtcau ljl997t) and 
iMathewson et all (|l992T ). iNavarro fc Steinmet j (|200dh 
calculated the specific angular momenta vs. rotational 
velocity for late-type spiral galaxies and compared them 
to numerical simulations. The discs were assumed to 
follow an exponential profile for which the peak rotational 
velocity, « ro t,2.2, occurs at r = 2.2rd, and it follows that 
jbar = 2rdWrot,2.2- This assumption can be misleading 
when comparing to simulated galaxies as the true v rot - 
peak can be significantly underestimated in th e case of 
bulge- dominated disc galaxies. The sample of ICourteaul 
(119971 ) concerned Sb-Sc galaxies for which B/D is low 
and a dominating exponential disc assumption is roughly 
valid. The difference in measured and estimated angular 
momentum content makes it difficult t o compare simulated 
and observed galaxies, as discuss ed in lAbadi et al.l (2003a) 
and IPiontek fc Steinmetzl (|2009bT ) . A simulated galaxy can 
be considered as a successful realizations of a late-type 
(Sb-Sc/Sd) galaxy if the estimated and measured angular 
momenta are in agreement. 

Focusing on the n = 0.1 cm -3 suite, we find that the 
nOlel and nOlelML simulations are in good agreement 
with the observed galaxies, both when analysed using Eg . 1111 
and the exponential disc approximation. Typical measured 
and estimated values are here jbar ~ 2000 km s _1 kpc and 
jbar ~ 2750kms~ 1 kpc respectively. As eg is increased, the 
calculated angular momentum decreases. The ee = 2 per 
cent simulations are still a part of the observed scatter while 
higher values create more significant outliers in the observed 
distribution. When using the exponential disc approxima- 
tion, all simulated galaxies are in good agreement with the 
observed data as the velocities are quite comparable at larger 
radii, and for the fact that the discs, although less extended, 
have larger rd in the higher efficiency cases (see Fig. [5}. We 
conclude that an angular momentum reservoir comparable 
to Sb/Sc galaxies have been reproduced for the baryons in 
the case of low SFE (i.e. eg = 1 per cent). 

The lack of correlation between Esnu and the bary- 
onic angular momentum content might come as a surprise. 
However, while the B/D ratio decreases for large super- 
nova energy injections, the actual disc mass changes lit- 
tle, and is ~ 6 x 10 10 M© for all eg — 5 per cent simu- 
lations. As the net contribution of the bulge to the angu- 
lar momentum content is roughly zero, similar jbar is to 
be expected. All eg = 5 per cent simulations have measured 
jbar ~ 1300—1450 kms _1 kpc which is close to the estimated 
jbar ~ 1600kms _1 kpc in n01e5SN5. 

In summary, the largest measured baryonic specific an- 
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Figur e 12. The i-ba nd Tully-Fisher relationship from the SDSS 
[|Pizagno et al.ll2007h . We show the observed average (solid line), 
lcr (dashed line) and 2<r (dotted line) relation. The symbols are 
results from our simulated galaxies, which use no = 0.1 cm -3 and 
eg = 1 (black symbols), 2 (red symbols) and 5 (blue symbols) per 
cent. 



gular momentum reservoir can be found in simulations using 
eg = 1 per cent due to a massive disc component, regard- 
less of including feedback or not. At higher efficiencies, jhar 
decreases, again regardless of feedback. 



7.2 The Tully-Fisher relationship 



The p hotometric Tully-Fisher (TF) relation (|Tullv fc Fisherl 
1 19771 ) links the characteristic rotational velocity of a galaxy 
with its total absolute magnitude. This correlation holds in 
all typical photo metric bands but wi th variation in func- 
tional form (e.g. Pizag no e t al. 2007) . Early atte mpts in 
forming realistic galaxies (e.g. lAbadi et al] l2003al ) showed 
off-sets in the observed relation owing to the forma- 
tion of very concentrated bulge-dominated galaxies with a 
low star formation activity at late times. Their velocity- 
magnitude relation h ad more in common with SO galaxies 
ijMathieu et alJl2002h . Recent work seems to have improved 
on these results by SN feedback regulated star fo rmation 
IIGovernato et al ] |2007l ; |Piontek fc Steinmet3l2009bl ). These 
studies place galaxies closer to the observed relation, but this 
is in part achieved by circumventing the large measured v TO t 
(caused by the dominant bulge) by using the exponential 
disc assumption discussed above (but see iGovernato et al] 
2009), i.e. using w ro t,2.2. Observationally, the measured quan- 
tity is often half of the HI velocity width at 20(VK2o) or 
50 (VK50) per cent of the peak intensity. 

In Fig.[l2]we present the measured i-band magnitudes 
of several of our simulated galaxies as a function of their 
peak rotational velocities measured from the gas compo- 
nent. Thes e are compared to t he observed TF relation from 
the SDSS jPizagno et alj|2007f ). Pizagno et al. measured the 
velocity at a radius containing 80 per cent of the i-band 
flux. This measure (Vso) is equivalent to measuring v TO t at 
~ 3rd for a pure exponential disc. By using the true peak 



of Vrot, we provide an absolute lower limit to the agreement 
with observations, and can clearly separate disc and bulge- 
dominated galaxies. We note that the low efficiency models 
agree well with the average data, regardless of adopted feed- 
back scheme, and even without. At higher efficiencies, the 
discs are off-set by more than 2a, mostly due to their peaked 
rotation curves. In these circumstances, the inclusion of ad- 
ditional recycling via SNIa and stellar mass-loss increases 
the magnitudes by ~ 0.5 dex in n01e2ML and n01e5ML 
simulations. The n01e5NFB simulation is brighter than the 
corresponding simulations including feedback due to exclu- 
sion of metal enrichment, leading to less efficient cooling and 
more gas left to form stars at later times. Allowing for en- 
richment without any energy deposition demonstrates this 
fact (see figure). From a photometric TF point of view, the 
eg = 5 per cent discs corresp ond to SO systems or early type 
spirals |Mathieu et al]|2002l ). 

As described in the previous sections, v rot and E s in the 
n01e5SN5 simulation agrees fairly well with the disc values 
found in nOlel. The strong feedback brings the galaxy closer 
to the observed values but the absolute magnitude is still 
lower than in the eg = 1 per cent simulations. The SFH in 
Fig.[8]tells us why: after z = 1 the SFR is lower in n01e5SN5 
compared to nOlel by almost a factor of 2 (even though the 
2 = values agree) due to strong gas expulsion, resulting in 
a less bright disc by ~ 1/3 dex in i-band magnitude. 

As for the specific angular momentum analysis, adopt- 
ing the "rot, 2. 2 measure (or Vso), all discs would agree sta- 
tistically with the observed TF relation, especially when in- 
cluding SNIa feedback and stellar mass-loss. 

Similar to the photometric T F is the 'baryonic TF' rela- 
tion (jMcGaueh et al.ll2000l . [20T0I ) which links characteristic 
rotation velocity with total galaxy baryonic mass. The bary- 
onic TF relation therefor accounts for the fact that less mas- 
sive galaxies are more gas rich, and their stars only account 
for a small fraction of the total disc mass. The same conclu- 
sion as above holds for the baryonic TF: the low-efficiency 
simulations agrees well with the observations. As the bary- 
onic masses for the discs are not strongly affected, even in 
the case of extreme feedback (n01e5SN5) the data points 
shift only with the increase of the v TO t peak. As for the pho- 
tometric TF, using "rot, 2. 2 puts all galaxies on the observed 
relation. 



7.3 The EsFR-Egas relation 

The most famous study of the globally averaged relationship 
betw een the star form ation rate and gas surface density is 
from lKennicuttl l| 19981 ) (from now on K98), where a sample 
of 61 nearby normal spiral galaxies and 36 infrared-selected 
starburst galaxies were considered. Assuming a Schmidt-law 
of the form 



Esfr = a 



1M pc- 



(12) 



the full sample yielded a = (2.5 ± 0.7) x 10 4 and 
N = 1.40 ± 0.15. The SsFR-S gas relation has been studied 



by many authors (e.g. Wong fc Blitz 20021: [Misiriotis et al] 
l2006l: iKennicutt et al]|2007l : Schuster et al ll2007D . both lo- 
cally and globally, using different star formation tracers 
and galaxy samples. A large range of power-law indices 
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Figure 13. Ssfr vs. S g for the resulting discs using eff = 1 (nOlel), 2 (n01e2) and 5 (n01e5) per cent at z = (left-hand p anel) and at 
z = 3 (right-hand panel). The filled circles are radial data for 7 spiral galaxies from the THINGS survey llBigiel et al.ll20o3) . where S gas 
includes the contribution from helium (S gas = 1.36Shi+H9)- The data points represent, from lightest to darkest, > 1, > 5, > 10, > 20 
and > 30 detections. The vertical dotted lines separate regions where different star formation laws are conjectured to apply (see text). 
Diagonal dotted lines show lines of constant SFE=SsFR/S gas , indicating the l evel of Ssfr need ed to consume 1, 10 and 100 per cent of 
the gas reservoir in 10 8 ye ars. The solid black li ne is the average relation from iKennicuttJ l|l998l) . The z = 3 observations approximately 
populate the region of the lWolfe fc Chenl i2006h observations. THINGS data courtesy of F. Bigiel. 
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Figure 14. Ssfr. vs. S g for the resulting discs using a high star formation efficiency (eff = 5 per cent), but with different SN feedback 
strengths: E S nii = 10 51 erg (n01e5), 2 X 10 51 erg (n01e5SN2), 5 X 10 51 erg (n01e5SN5) as well as with zero SNII feedback energy but 
metal enrichment (n01e5NFBmet). The panels show the results at z = (left) and at z = 3 (right). The lines and symbols are described 
in the caption of Fig. 1131 SN feedback has little effect on the SsFR-S g relation at z = but does affect the high redshift relation, although 
only for very large energy injections (i?SNll > 2 X 10 erg). 
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(N as 1 — 3) have been found, suggesting that either different 
SF laws exist in different galaxies or that N is very sensitive 
to systematic differences in methodology. IBigiel et all (2008) 
presented a comprehensive analysis of the EsFR-E gas rela- 
tionship using multifrequency data of 7 spiral galaxies and 
f f late-type and dwarf galaxies. The analysis pointed to a 
great variation within the sample and a markedly different 
functional behaviour in atomic-and molecular-dominated 
gas. 

The THINGS data of Bigiel et al., relevant for spirals, 
as well as the K98 law fEq . 112(1 . is reproduced in the left- 
hand panel in Fig.[l3]together with the azimuthally averaged 
(Ar = 540 pc) data from nOlel, n01e2 and n01e5. For the 
calculation of Esfr we only consider stars younger than 50 
Myr. At a given value of E gaa we find a clear trend of higher 
Esfr values for higher eg. All simulations are compatible 
with the range of observed values, having the same func- 
tional behaviour but with an off-set. We note that only the 
disc in the n01e5 simulations is compatible with the K98 
relation. The nOlel simulation is on the low side but can 
still statistically be associated with one of the THINGS spi- 
ral galaxies. However, at high redshift the argument can be 
reversed, as can be seen in the right- hand panel of F ig.fTiH 
The observations of DLAs at z ~ 3 by I Wolfe fc Chenl (2006) 
are typically an order of magnitude lower than the K98 rela- 
tion, agreeing only with measurement of the low density en- 
vironment of the discs in our low efficiency simulations. This 
trend is also pre dicted by simulations includ ing treatment 
of H2 formation (|Gnedin fc Kravtsov|[2010bl lah. In essence, 
while nOlel is on the low side at z = 0, it is consistent 
with high redshift observations and the reverse argument is 
valid for n01e5. A higher efficiency is acceptable at lower 
redshift, and is predicted due to e.g. higher gas metallicity. 
As the bulge component is assembled at high redshift, the 
efficiency of star formation during this epoch is crucial in 
setting the morphology of the galaxy. 

The same analysis is performed for the feedback test 
suite in Section|6] and shown in Fig. 1141 At z = 0, all sim- 
ulations show a similar functional behaviour, but with a 
weak trend of lower Esfr as -Esnii is increased, while re- 
maining comparable to the K98 law. At z = 3, a slightly 
greater effect is found, but only for very large energy injec- 
tions (-Esnii > 2 x 10 5 erg). The extreme case of Esnii = 
5 x 10 51 erg (n01e5SN5) is comparable to a lowering star 
formation efficiency to eg = 2 per cent (n01e2). None of the 
strong feedback simulations regulate star formation enough 
to reproduce the low Esfr values found for eg = 1 per cent 
(nOlel). This z ~ 3 ins ensitivity of the K -S relation to feed- 
back was also found bv lKravtsovl (|2003h . 



8 DISCUSSION AND CONCLUSIONS 

In this paper we have presented a set of AMR simulations 
studying the assembly of large Milky Way-like disc galaxies. 
The self-consistent formation of a late-type disc galaxy has 
remained elusive in the field of numerical galaxy formation, 
mainly due to the strong loss of angular momentum in the 
galaxy assembly process. A popular solution to this problem 
is to regulate star formation at high redshift via supernova 
explosions that drive galactic winds, transporting material 



out of star-forming regions hence lowering the local star for- 
mation rate. 

We have investigated the plausibility of this mechanism 
in comparison to a small scale (~ 100 pc) physical approach 
where star formation is made inefficient by modifying the 
Schmidt-law star formation normalization. In a very crude 
way, this mimics unresolved physics such as H2 formation, 
small scale turbulence and radiative effects. We find that the 
Schmidt-law efficiency of star formation is far more success- 
ful way of regulating star formation towards realistic galax- 
ies than what can be achieved via supernova feedback. Our 
most successful models reproduce Milky Way galaxies with 
flat rotation curves, where the small bulge component is 
formed via secular processes. The main conclusions of this 
work can be summarized as follows. 

(i) Disk characteristics such as E*(r), E gas (r), v TO t{r) and 
B/D strongly depend on the choice of star formation effi- 
ciency per free-fall time, eg. The parameter will essentially 
set the mode of global star formation, hence governing the 
final spiral Hubble type, where low efficiencies of eg ~ 1 
per cent render discs of Sb or Sbc type, while eg — 5 
per cent moves the discs closer to Sa/SO types. Simulations 
at low efficiencies agre e well with observational constraint s 
on disc characteristics (|Courteaulll9"97l ; iGnedin et af1l2007f ) , 
as well as the angular mom entum content of disc galaxies 
I Nava rro & Stcinmctz 20Q0|) , the Tully- Fisher r elationship 
( Pizagno et al.ll2007l) an d the EsFR-E gas relation (|Kennicuttl 
Il998t IBigiel et al.ll200sl ). The origin of the successful Milky 
Way-like galaxy formation is a well motivated suppression 
of star formation at z ~ 3, the epoch at which the violent 
assembly process would form a slowly rotating bulge in case 
of efficient star formation. 

(ii) Supernova feedback does not regulate star formation 
efficiently at low input energies. Only when the injected en- 
ergy per supernova event is five times the canonical value, i.e. 
5xl0 51 erg, do we find lower and more realistic B/D ratios in 
the simulations tuned to the standard lKennicuttl l|l998l ) star 
formation law, leading to a flatter rotational velocity profile, 
hence resembling the galaxies formed without strong feed- 
back but with a low Schmidt-law efficiency. This comes at 
the cost of a significantly distorted gas disc at 2 = 0, as 
well as a less bright stellar disc as gas is expelled into the 
halo, leaving less fuel for star formation at late times. In 
essence, we find that changes in eg can play a much greater 
role in shaping a spiral galaxy than gas redistribution via 
supernovae-driven winds. 

It is plausible that at very high resolution, or using a dras- 
tically different recipe of supernovae feedback, lower values 
of -Esnii may be successful in regulating the SFE. If so, it 
will still need to mimic the low efficiency on scales of a few 
100 pc which, as argued in this work, can be absorbed by 
the en-term. 

(iii) If the star formation efficiency paramete r is tuned 
to m atch the standard z = K-S data ( Ken nicuttl 
19981) . i.e. requirin g on the order of eg > 5 percent (e.g. 
Stinson et al.l 120061 ) . star formation is likely to be overes- 



timated at high redshift (z = 3) whe re the amplitudes of 
Esfr are an order of magn itude lower (|Wolfe fc Chen|[2C306l ; 
IGnedin fc Kravtsov|[2010bl ). All efficiencies studied in this 
work (eg = 1 — 5 per cent) are compatible wit h modern 
data of the THINGS survey (|Bigiel et al.l [2008) but only 
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when eg ~ 1 per cent can the constraints from z = 3 data 
be met and late-type, disc dominated systems form. As the 
true SFE varies in space and time, being depe ndent on small 
scale physics governing H2 formation (see e.g. iGnedin et al.l 
2009), present day simulations based on single valued effi- 
ciency parameter have little predictive power. 



We argue (see also IGnedin et aTJ [2009) that the re- 
sults presented in this paper indicate that other processes 
in the ISM in addition to, or in conjunction with, super- 
nova feedback are important in explaining the evolution of 
the galaxy population, as well as regulating observed disc 
sizes. Some form of outflow process must be res ponsible for 
enriching the IGM ([Oppenheimer fc Pavel l2006h . which to- 
gether with an inefficient star forma tion might explain the 
faint end of the stellar m ass function (|Somerville fc Primackl 
Il999l : iKeres et al . 2009). The same argument can be used 
for the mass-metallic i ty rel ationship ([Brooks et al] 12007? ) . 
although pTassis et all l[200S?) demonstrated that it could be 
reproduced without supernova-driven outflows. Galaxies of 
masses considered in this work are situated at the knee of 
the stellar mass function, where the o bserved and simul ated 
functions (even without feedback; see lKeres et aljl2009? ) are 
in closest agreement. This circumstance might explain why 
even our simulations without feedback resulted in realistic 
discs. At this galaxy mass, supernova driven winds can- 
not escape the deep potential well, and are impeded by 
the hot halo. On the other hand, AGN feedback, which re- 
cently has been introduced into galaxy formation simula- 
tions ibi Matteo et all [2005). is probably not relevant for 
the Milky Way since the black hole might not be massive 
enough for efficient AGN radio-heating. At higher masses, 
and/or at high redshift, the inclusion of AGN is probably 
necessary to correctly reproduce the observed abundances 
and stellar masses. This is the greatest uncertainty of our 
work, which we leave for a future study. 

The way in which galaxies p opulate dark matter h aloes 
is an important topic, see e.g. button et ail \20Vj ) and 
references within for a compilation of r ecent observationa l 
data an d theoretical work. Rece ntly. iGuo et al.l (12010? ) 
[see also iMoster et all (12010? ) and behroozi et"aiT < |2010h ] 
matched dark matter halo mass function from cosmolog- 
ical N— body simulations to the stellar mass function of 
the galaxies from the SDSS ([Li fc White! 12009? ) . This anal- 
ysis yields the required galaxy formation efficiency, r\ = 
(M*/Mhaio)(fim/fib), i-e. what fraction of the universal 
baryons that must have condensed into stars at a given halo 
mass. In our "best-case" model (nOlelML, see Table[T|, the 
total stellar and dark matter halo virial mass is ~ 10 11 M 
and ~ 10 12 Mq respectively. This results in a stellar frac- 
tion of 10 per cent, which corresponds to almost 60 per cent 
of the cosmic baryon fraction. The rest of the baryons re- 
side in the stellar halo, gaseous disc and ionized gas halo. At 
this halo mass, abundance matching requires that the stellar 
disc accounts for only ~ 20 per cent of the cosmic baryon 
fraction, i.e. a factor of three lower. Similar discrepancies 
exist in all modern work of numerical galaxy formation 
(lAbadi et al]|2003al; lOkamoto etafl 120051; [Governato et all 



12007? ; IScannapieco et al]|2009l ; |Piontek fc Steinmetzll2009b? l. 
and its origin is not yet know, although AGN is a compelling 
mechanism at the high mass end, as discussed above. This 
issue is the topic of a follow-up paper in preparation. 



behroozi et al] (?2010T ) performed a comprehensive anal- 
ysis of abundance matching, accounting for systematic errors 
in e.g. the stellar mass estimates, the halo mass function, 
cosmology etc. Our simulated galaxy formation efficiencies 
would be in ~ 2a agreement with their result (see their Fig. 
11). We note that our own Galaxy and M31 also might be 
strong outliers in thi s analysis, considering the inferred r \ 
from mass modelling |Klypin et al.ll2002l ; ISeigar et alj|2008?) 
as we ll as via recent MW halo mass estimates |Xue et al] 
12008? ). 

Abundance matching is insensitive to the actual Hub- 
ble types of the galaxies, forcing all galaxies of a specific 
mass to be linked to only one halo mass. At the stellar 
mass scale of the Milky Way (5 - 7 x IP 1 " M B ), only ~ 25 
per c ent of galaxies are of Sb/Sbc type l|Nair fc Abraham! 
12010? ). It is plausible that the more active merger histories 
associated with ellipticals and early type disc galaxies have 
led to a stronger mass expulsion, via e.g. AGN, in com- 
parison to the more disc dominated counterparts. In this 
scenario, late-type discs are expected to be outliers in the 
galaxy formation efficiency vs. stellar mass relation, consid- 
ering the strong bias towards early type systems. A detailed 
sub-division into Hubble types has not yet been performed 
when matching galaxies to haloes, although color separa- 
tions into red and b lue systems have been m ade in studies 
using weak- lensing ([Mandelbaum et alj|2006? ) and satellite 
kinematics ([More et al.ll2010? ). These studies indicate a dif- 
ferent galaxy formation efficiency for galaxies similar in mass 
to the Milky Way; a late-type galaxy is associated with a 
halo of ~ 0.5 dex lower halo mass c ompared to an eq ually 
massive early type (see e.g. fig. 11 in I More et al.l[201ol ). 

Understanding, from a numerical perspective, the 
spread of baryon fractions across dark matter haloes of dif- 
ferent masses, accretion histories and environments is a com- 
plicated problem, and will require a large sample of high- 
resolution simulations, which we leave for a future investi- 
gation. 
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APPENDIX A: EFFECT OF RESOLUTION 

The numerical setup in the SR5-nlelML simulation is iden- 
tical to SR6-nlelML, apart from having twice as high spa- 
tial resolution and hence eight times the mass resolution (see 
Table[l|. In Fig. lAll we plot the resulting stellar surface den- 
sities and circular velocities from the simulations. We find 
very little difference apart from more gas being consumed 
in SR5-nlelML resulting in a more massive disc (by ~ 35 
per cent) at higher resolution and the bulge mass is lowered 
by ~ 20 per cent. The dark halo is less contracted but as 
M(< r) remains roughly the same, the circular velocities 
changes little. 

Why do we not see any striking differences in the higher 
resolution simulation? At the adopted resolutions we do not 
uncover more physics with only a factor of two increase 
in spatial resolution: we are still under-resolving the scale 
height and are not resolving individual star-forming clouds. 
At much higher resolution this conclusion will change and 
the numerical parameters will have a different meaning. 
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Figure Al. Stellar surface densities (top) and circular velocities 
(bottom) in SR6-nlelML (black line) and SR5-nlelML (red line). 



can then be calculated as 
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The true star-forming scale height of cold molecular gas is 
on the order of ~ 10 pc, and for the neutral atomic ISM 
it is ~ 100 pc. Galaxy formation simulations under-resolve 
this structure, and it is safe to assume that the effective 
numerical scale height h ~ Ax, until Aa; goes below a few 
tens of parsecs. Under this assumption the integration limits 
in the sech 2 term above will always be from — z/h = —0.5 
to z/h — 0.5. Eq. 5 can now be simplified to 
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which in units of particles per cm 3 this is approximately 
f Md \ 3exp(-r/r d ) _ 3 
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Axr\ 
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We can rearrange this equation and write down the following 
useful relationship: 



r — rd In 
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(B7) 



This equation tells us at which radius, given a cell size, we 
cross the density uax- In Section lSTTI we argued that using 
n — 1 cm -3 and Aa; = 340 pc resulted in less realistic galax- 
ies as star formation was "missing" at large radii. Inserting 
these values into Eq. lB7l gives us r ~ 6—8 kpc for any reason- 
able choice of r^, confirming the arguments of Section [5~T1 
Adopting n = 0.1 cm" 
a factor of 2.3. 



increases the star-forming radius by 



APPENDIX B: RESOLVING STAR 
FORMATION 

In this appendix we derive a simple relationship of the re- 
quired resolution (Ax) and star formation density threshold 
(no) required to capture star formation in an extended disc 
at the current epoch. Let us assume that a gas disc at z = 
follows an exponential radial density profile with a sech 2 
vertical profile, i.e. 

p(r,z) = p sech 2 (-z/h) exp(-r/r d ). (Bl) 

The mass of this profile can be integrated from Eq. lBll to 
give us the characteristic scale disc density 

po = M d /4nhr 2 d . (B2) 

Assume that this density distribution is coarsened on a reg- 
ular mesh with a cell size Ax, and that the disc is aligned to 
a mesh axis. The physical density in a central strip of cells 



